  module sparse
! BiCGStab solver
  use prec

  contains

  subroutine bicgstab(A,b,n,x0,x,maxit,tol,relres,niter)

  IMPLICIT NONE

!  INTEGER, PARAMETER :: dp = KIND(1.0D0)
  INTEGER :: n
  INTEGER :: niter,maxit
  REAL(dp) :: tol
  REAl(dp), DIMENSION(n,n) :: A
  REAl(dp), DIMENSION(n) :: b
  REAl(dp), DIMENSION(n) :: x0
  REAl(dp), DIMENSION(n) :: x

!
! Locals
!
  INTEGER :: i,j
  REAL(dp) :: rho, alpha, om, beta, om_old, rho_old, suma, omim, rhoim, normr, normb, relres
  REAL(dp),DIMENSION(n) :: r, r0h, v, p, s, t, r_old, v_old, p_old, x_old
! Initialize
  do i = 1,n
  r(i) = b(i)-sum(A(i,:)*x0(:)) ! vectors
  enddo
  r0h = r
  rho = 1. ;alpha = 1.;  om = 1.0_dp
  v = 1.; p = 0.0_dp ! vectors
  normb = sqrt(sum(b*b))

! Main loop
  do j = 1,maxit
  niter = j

! Shift - save values from previous iteration
  r_old = r
  v_old = v
  p_old = p
  x_old = x
  om_old = om
  rho_old = rho

  rho = sum(r0h*r_old)
  rhoim = 1./rho_old
  omim = 1./om_old
  beta = rho*rhoim*alpha*omim
  p = r_old + beta*(p_old - om_old*v_old)
  do i = 1,n
  v(i) = sum(A(i,:)*p(:))
  enddo
  suma = 1./sum(r0h*v)
  alpha = rho*suma
  s = r_old - alpha*v
  do i = 1,n
  t(i) = sum(A(i,:)*s(:))
  enddo
  suma = 1./sum(t*t)
  om = sum(t*s)*suma
  x = x_old + alpha*p + om*s 
  r = s - om*t
  normr = sqrt(sum(r*r))
  relres = normr/normb
  if(relres < tol.or.j>=maxit) exit
  enddo

  return
  end subroutine

      function solve(a,b,n) result(u)
!     driver for bicgstab routine
!      use prec
!      use sparse
      integer            :: niter = 1000
      real(dp) :: tol = 1e-6
      real(dp) :: relres
      real(dp), dimension(n,n), intent(in) :: A
      real(dp), dimension(n), intent(in) :: b
      real(dp), dimension(n) :: u0 ! initial solution vec
      real(dp), dimension(n) :: u  ! solution vec
      u = 0.0_dp
      u0 = 0.0_dp

      call bicgstab(A,b,n,u0,u,1000,tol,relres,niter)

      end function solve

  end module sparse
